Thesis - Model Building 1
Jacob Levine
2/6/2018
## load packages
library(mgcv)
library(plyr)
library(visreg)
library(ggplot2)
library(MuMIn)
library(reticulate)Model Building
## read in the main data
allom_data <- read.csv('~/Desktop/References_allometry/WMNF_allom5.csv', header = TRUE)
str(allom_data)## 'data.frame': 318 obs. of 21 variables:
## $ SPP4 : Factor w/ 8 levels "ACPE","ACRU",..: 2 2 2 3 3 3 4 4 4 5 ...
## $ ID : Factor w/ 68 levels "1","10","11",..: 55 55 55 55 55 55 55 55 55 55 ...
## $ ELEVB : int NA NA NA NA NA NA NA NA NA NA ...
## $ TREEAGE : int NA NA NA NA NA NA NA NA NA NA ...
## $ STANDAGE : int 14 14 14 14 14 14 14 14 14 14 ...
## $ LOCATION : Factor w/ 3 levels "BEF","HBEF","WMNF": 1 1 1 1 1 1 1 1 1 1 ...
## $ DataSetID : Factor w/ 6 levels "BattlesGB","BattlesRock",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ SOURCE : Factor w/ 5 levels "BattlesJ","FaheyT",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ DBHcm : num 2.2 3.4 4.4 2 2.9 3.9 2.4 3 3.8 3.3 ...
## $ HTcm : num 554 520 685 466 647 569 455 430 625 840 ...
## $ PVcm3 : num 1053 2361 5208 732 2137 ...
## $ LEAFg : int 69 179 415 89 279 341 97 280 330 142 ...
## $ TWIGg : int NA NA NA NA NA NA NA NA NA NA ...
## $ CTLGRg : num NA NA NA NA NA NA NA NA NA NA ...
## $ LIVEBRANCHg: num 281.5 376.8 779.9 97.4 311.6 ...
## $ BOLEBARKg : num NA NA NA NA NA NA NA NA NA NA ...
## $ BOLEWOODg : num NA NA NA NA NA NA NA NA NA NA ...
## $ BOLETOTg : num NA NA NA NA NA NA NA NA NA NA ...
## $ ABOVEg : num 848 1986 4283 745 1916 ...
## $ DEADWOODg : num NA NA NA NA NA NA NA NA NA NA ...
## $ TOTROOTg : num NA NA NA NA NA NA NA NA NA NA ...allom_data[, "STANDAGE"] <- as.factor(allom_data[, "STANDAGE"])## Remove the data from the 98 year old stand so that all of the data are comparable:
allom_data_Y <- subset(allom_data, STANDAGE != "98years")
str(allom_data_Y)## 'data.frame': 318 obs. of 21 variables:
## $ SPP4 : Factor w/ 8 levels "ACPE","ACRU",..: 2 2 2 3 3 3 4 4 4 5 ...
## $ ID : Factor w/ 68 levels "1","10","11",..: 55 55 55 55 55 55 55 55 55 55 ...
## $ ELEVB : int NA NA NA NA NA NA NA NA NA NA ...
## $ TREEAGE : int NA NA NA NA NA NA NA NA NA NA ...
## $ STANDAGE : Factor w/ 7 levels "14","16","17",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ LOCATION : Factor w/ 3 levels "BEF","HBEF","WMNF": 1 1 1 1 1 1 1 1 1 1 ...
## $ DataSetID : Factor w/ 6 levels "BattlesGB","BattlesRock",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ SOURCE : Factor w/ 5 levels "BattlesJ","FaheyT",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ DBHcm : num 2.2 3.4 4.4 2 2.9 3.9 2.4 3 3.8 3.3 ...
## $ HTcm : num 554 520 685 466 647 569 455 430 625 840 ...
## $ PVcm3 : num 1053 2361 5208 732 2137 ...
## $ LEAFg : int 69 179 415 89 279 341 97 280 330 142 ...
## $ TWIGg : int NA NA NA NA NA NA NA NA NA NA ...
## $ CTLGRg : num NA NA NA NA NA NA NA NA NA NA ...
## $ LIVEBRANCHg: num 281.5 376.8 779.9 97.4 311.6 ...
## $ BOLEBARKg : num NA NA NA NA NA NA NA NA NA NA ...
## $ BOLEWOODg : num NA NA NA NA NA NA NA NA NA NA ...
## $ BOLETOTg : num NA NA NA NA NA NA NA NA NA NA ...
## $ ABOVEg : num 848 1986 4283 745 1916 ...
## $ DEADWOODg : num NA NA NA NA NA NA NA NA NA NA ...
## $ TOTROOTg : num NA NA NA NA NA NA NA NA NA NA ...max(allom_data_Y$DBHcm)## [1] 66## get subset of data which includes trees in 98year stand which are smaller than the ## largest trees in the younger stands.
allom_data_trunc <- subset(allom_data, DBHcm <= 66)## log transform aboveground biomass:
allom_data_trunc$logABOVEg <- log(allom_data_trunc$ABOVEg)
ggplot(allom_data_trunc, aes(x = DBHcm, y = logABOVEg)) +
geom_jitter() +
geom_smooth(size = 0.5)
## lets try taking the log transformation of both variables:
allom_data_trunc$logDBHcm <- log(allom_data_trunc$DBHcm)
ggplot(allom_data_trunc, aes(x = logDBHcm, y = logABOVEg)) +
geom_jitter() +
geom_smooth(size = 0.5)
## Create SPPSTAGE in order to fully cross in GAM
for (i in 1:length(allom_data_trunc$SPP4)) {
allom_data_trunc$SPPSTAGE[i] <- paste(allom_data_trunc$SPP4[i], allom_data_trunc$STANDAGE[i])
}
allom_data_trunc$SPPSTAGE <- as.factor(allom_data_trunc$SPPSTAGE)
## log transform HTcm:
allom_data_trunc$logHTcm <- log(allom_data_trunc$HTcm)
## log transform PVcm3
allom_data_trunc$logPVcm3 <- log(allom_data_trunc$PVcm3)
## remove NA values of biomass:
allom_data_trunc <- subset(allom_data_trunc, logABOVEg != "NA")## look at how many data points we have for each species. ACRU and PIRU seem low.
## Let's create a dataset with them removed to compare to one with them included.
count(allom_data_trunc$SPP4)## x freq
## 1 ACPE 40
## 2 ACRU 12
## 3 ACSA 49
## 4 BEAL 62
## 5 BEPA 41
## 6 FAGR 48
## 7 PIRU 15
## 8 PRPE 49allom_data_truncSPP4 <- subset(allom_data_trunc, SPP4 != "ACRU")
allom_data_truncSPP4 <- subset(allom_data_truncSPP4, SPP4 != "PIRU")
count(allom_data_truncSPP4$SPP4)## x freq
## 1 ACPE 40
## 2 ACSA 49
## 3 BEAL 62
## 4 BEPA 41
## 5 FAGR 48
## 6 PRPE 49## make another subset of the data with only SPPSTAGE levels that have more than 10 data points.
SPPSTAGE_counts <- data.frame(count(allom_data_trunc$SPPSTAGE))
SPPSTAGE_counts <- subset(SPPSTAGE_counts, freq >= 10)
SPPSTAGE_list <- as.factor(SPPSTAGE_counts$x)
allom_data_truncSPPSTAGE <- subset(allom_data_trunc, (allom_data_trunc$SPPSTAGE %in% SPPSTAGE_list) == TRUE)
allom_data_truncSPPSTAGE$SPPSTAGE## [1] BEAL 17 BEAL 17 BEAL 17 BEAL 17 BEAL 17 BEAL 17 BEPA 17 BEPA 17
## [9] BEPA 17 BEPA 17 BEPA 17 BEPA 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17
## [17] PRPE 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17
## [25] PRPE 17 PRPE 17 PRPE 17 ACPE 17 ACPE 17 ACPE 17 ACPE 17 ACPE 17
## [33] ACPE 17 ACPE 17 ACPE 17 ACPE 17 ACPE 17 ACSA 17 ACSA 17 ACSA 17
## [41] ACSA 17 ACSA 17 ACSA 17 ACSA 17 ACSA 17 ACSA 17 ACSA 17 ACSA 17
## [49] BEAL 17 BEAL 17 BEAL 17 BEAL 17 BEAL 17 BEAL 17 BEAL 17 BEAL 17
## [57] BEAL 17 BEAL 17 BEAL 17 BEPA 17 BEPA 17 BEPA 17 BEPA 17 BEPA 17
## [65] BEPA 17 BEPA 17 BEPA 17 BEPA 17 BEPA 17 BEPA 17 FAGR 17 FAGR 17
## [73] FAGR 17 FAGR 17 FAGR 17 FAGR 17 FAGR 17 FAGR 17 FAGR 17 FAGR 17
## [81] FAGR 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17
## [89] PRPE 17 PRPE 17 PRPE 17 ACPE 26 ACPE 26 ACPE 26 ACPE 26 ACPE 26
## [97] ACPE 26 ACPE 26 ACPE 26 ACPE 26 ACPE 26 ACPE 26 ACPE 26 ACPE 26
## [105] ACPE 26 ACPE 26 BEAL 26 BEAL 26 BEAL 26 BEAL 26 BEAL 26 BEAL 26
## [113] BEAL 26 BEAL 26 BEPA 26 BEPA 26 BEPA 26 BEPA 26 BEPA 26 BEPA 26
## [121] BEPA 26 BEPA 26 PRPE 26 PRPE 26 PRPE 26 PRPE 26 PRPE 26 PRPE 26
## [129] PRPE 26 PRPE 26 BEAL 26 BEAL 26 BEAL 26 BEPA 26 BEPA 26 BEPA 26
## [137] PRPE 26 PRPE 26 PRPE 26 ACPE 98 ACPE 98 ACPE 98 ACPE 98 ACPE 98
## [145] ACPE 98 ACPE 98 ACPE 98 ACPE 98 ACPE 98 ACPE 98 ACPE 98 ACPE 98
## [153] ACPE 98 ACPE 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98
## [161] ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98
## [169] ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 BEAL 98
## [177] BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98
## [185] BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98
## [193] BEAL 98 BEAL 98 BEAL 98 BEAL 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98
## [201] FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98
## [209] FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98
## [217] FAGR 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98
## [225] PIRU 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98
## 39 Levels: ACPE 17 ACPE 26 ACPE 98 ACRU 14 ACRU 16 ACRU 26 ... PRPE 29## Try STANDAGE as numerical
allom_data_truncSTANDAGENUM <- allom_data_trunc
allom_data_truncSTANDAGENUM$STANDAGE <- as.numeric(allom_data_trunc$STANDAGE)Perhaps in the future we could keep the small trees only from the 98 year data as there is some basis to comparison, although its perhaps still too large of a gap to justify any conclusions.
It looks like we get something roughly linear with the loglog transformation of DBH vs BIOMASS. Based on a visual inspection alone, it appears that the heteroscedasticity issue has been at least significantly reduced.
Define Function to Create List of Models (linear and GAM):
A note: A more experienced/clever programmer could probably do this in 1/4 the lines of code I used. I am no such programmer, so this is likely the most overcomplicated and inefficient method of doing this. Oh well, maybe some day…
## comparing all possible combinations of linear models for the dataset using AIC and
## AICc:
## knowledge of the dredge function may have been useful/done the work for me for the gam
## models
options(na.action = "na.fail")
model1 <- lm(logABOVEg ~ logDBHcm * logHTcm * SPP4 * STANDAGE, data = allom_data_trunc)
model_sel_linear_AIC <- dredge(model1, rank = AIC)
model_sel_linear_AICc <- dredge(model1, rank = AICc)
## Function to write out all possible model terms, to be combined in later funcitons:
build_model_terms <- function(x, by, op) {
s_terms <- character()
int_terms <- character()
## for loops to write all possible smooth terms:
## l is dummy variable specifying inclusion of id = 1:
for (l in 1:2) {
if (l == 1) {
## k is a dummy variable indicating whether or not theres an interaction term:
for (k in 1:2) {
if (k == 1) {
for (i in 1:length(x)) {
for (j in 1:length(by)) {
s_termtext <- paste0("s(", x[i], ", k = 4", ", by =", by[j], " ,id = 1", ")")
s_terms <- c(s_terms, s_termtext)
}
}
}
else if (k == 2) {
for (j in 1:length(by)) {
s_termtext <- paste0("s(", x[1], ", ", x[2], ", k = 4", ", by =", by[j], " ,id = 1", ")")
s_terms <- c(s_terms, s_termtext)
}
}
}
}
else if (l == 2) {
## k is a dummy variable indicating whether or not theres an interaction term:
for (k in 1:2) {
if (k == 1) {
for (i in 1:length(x)) {
for (j in 1:length(by)) {
s_termtext <- paste0("s(", x[i], ", k = 4", ", by =", by[j], ")")
s_terms <- c(s_terms, s_termtext)
}
}
}
else if (k == 2) {
for (j in 1:length(by)) {
s_termtext <- paste0("s(", x[1], ", ", x[2], ", k = 4", ", by =", by[j], ")")
s_terms <- c(s_terms, s_termtext)
}
}
}
}
}
## k is a dummy variable indicating whether there are 1 or 2 intercept terms
for (k in 1:2) {
if (k == 1) {
for (i in 1:length(by)) {
by1 <- by[! by %in% " "]
if (!is.na(by1[i])) {
int_termtext <- paste0(by[i])
int_terms <- c(int_terms, int_termtext)
}
else {}
}
}
else if (k == 2) {
for (l in 1:2) {
int_termtext <- paste0(by[1], op[l], by[2])
int_terms <- c(int_terms, int_termtext)
}
}
}
model_terms <- list(int_terms, s_terms)
return(model_terms)
}
## Function builds models with one s() term - including those with two explanatory variables and an
## interaction:
build_models_1 <- function(y, x, by, int_terms, s_terms) {
models <- character()
for (i in 1:length(s_terms)) {
for (j in 1:length(int_terms)) {
## k is a dummy variable that indicates whether there are 1 or two intercept terms
for (k in 1:2) {
if (k == 1) {
## When we have only 1 intercept term:
if (!grepl("\\+", int_terms[j]) && !grepl("\\*", int_terms[j])) {
if (grepl(int_terms[j], s_terms[i]) && !grepl("SPPSTAGE", s_terms[i])) {
newmodeltext <- paste0(y, "~ ", s_terms[i], " + ", int_terms[j])
newmodel <- as.formula(parse(text = newmodeltext)[[1]])
models <- c(models, newmodel)
}
else if (grepl("SPPSTAGE", s_terms[i])) {
newmodeltext <- paste0(y, " ~ ", s_terms[i], " + ", "SPPSTAGE")
newmodel <- as.formula(parse(text = newmodeltext)[[1]])
models <- c(models, newmodel)
}
else if (grepl("by = )", s_terms[i]) && !grepl("SPPSTAGE", int_terms[j])) {
newmodeltext <- paste0(y, " ~ ", s_terms[i], " + ", int_terms[j])
newmodel <- as.formula(parse(text = newmodeltext)[[1]])
models <- c(models, newmodel)
}
else {}
}
else {}
}
else if (k == 2) {
if (grepl("SPPSTAGE", s_terms[i])) {}
else if (grepl("\\+", int_terms[j]) | grepl("\\*", int_terms[j])) {
newmodeltext <- paste0(y, " ~ ", s_terms[i], " + ", int_terms[j])
newmodel <- as.formula(parse(text = newmodeltext)[[1]])
models <- c(models, newmodel)
}
}
}
}
}
models <- unique(models)
return(models)
}
## Function to build models for case with two s() terms, no interactions:
## I got a bit carried away with this one I think.
build_models_2 <- function(y, x, by, int_terms, s_terms) {
models <- character()
## create a vector of DBH s() terms and one of HT s() terms
dbh_terms <- character()
ht_terms <- character()
for (i in 1:length(s_terms)) {
if (grepl("logDBHcm", s_terms[i]) && !grepl("logHTcm", s_terms[i])) {
dbh_terms <- c(dbh_terms, s_terms[i])
}
else if (grepl("logHTcm", s_terms[i]) && !grepl("logDBHcm", s_terms[i])) {
ht_terms <- c(ht_terms, s_terms[i])
}
else {}
}
## Make sure there are no duplicates
dbh_terms <- unique(dbh_terms)
ht_terms <- unique(ht_terms)
for (i in 1:length(dbh_terms)) {
for (l in 1:length(ht_terms)) {
for (j in 1:length(int_terms)) {
for (k in 1:2) {
if (k == 1) {
## When we have only 1 intercept term (except with special case of SPP4STANDAGE):
if (!grepl("\\+", int_terms[j]) && !grepl("\\*", int_terms[j])) {
## for the case without SPP4STANDAGE, require that both dbh_term and ht_term are represented in int_term:
if ((grepl(int_terms[j], dbh_terms[i]) && grepl(int_terms[j], ht_terms[l])) | (grepl(int_terms[j], paste0(dbh_terms[i], ht_terms[l])) && grepl("by = )", paste0(dbh_terms[i], ht_terms[l]))) && !grepl("SPPSTAGE", paste0(dbh_terms[i], ht_terms[l]))) {
newmodeltext <- paste0(y, " ~ ", dbh_terms[i], " + ", ht_terms[l], " + ", int_terms[j])
newmodel <- as.formula(parse(text = newmodeltext)[[1]])
models <- c(models, newmodel)
}
## for the case with SPPSTAGE:
else if (grepl("SPPSTAGE", paste0(dbh_terms[i], ht_terms[i]))) {
## for the case with SPP4STANDAGE and different by term, require that diff by term is represented in int_term
if (grepl(int_terms[j], paste0(dbh_terms[i], ht_terms[l])) && int_terms[j] != "SPPSTAGE") {
newmodeltext <- paste0(y, " ~ ", dbh_terms[i], " + ", ht_terms[l], " + ", "SPPSTAGE", " + ", int_terms[j])
newmodel <- as.formula(parse(text = newmodeltext)[[1]])
models <- c(models, newmodel)
}
## for the case with only SPPSTAGE as a by term, only SPPSTAGE as an int_term
else if (grepl(int_terms[j], dbh_terms[i]) && grepl(int_terms[j], ht_terms[l])) {
newmodeltext <- paste0(y, " ~ ", dbh_terms[i], " + ", ht_terms[l], " + ", "SPPSTAGE")
newmodel <- as.formula(parse(text = newmodeltext)[[1]])
models <- c(models, newmodel)
}
else {}
}
## for the case with empty by terms
else if (grepl("by = )", dbh_terms[i]) && grepl("by = )", ht_terms[l]) && !grepl("SPPSTAGE", int_terms[j])) {
newmodeltext <- paste0(y, " ~ ", dbh_terms[i], " + ", ht_terms[l], " + ", int_terms[j])
newmodel <- as.formula(parse(text = newmodeltext)[[1]])
models <- c(models, newmodel)
}
else {}
}
else {}
}
else if (k == 2) {
if (grepl("SPPSTAGE", paste0(int_terms[j], dbh_terms[i], ht_terms[l]))) {}
else if (grepl("\\+", int_terms[j]) | grepl("\\*", int_terms[j])) {
newmodeltext <- paste0(y, " ~ ", dbh_terms[i], " + ", ht_terms[l], " + ", int_terms[j])
newmodel <- as.formula(parse(text = newmodeltext)[[1]])
models <- c(models, newmodel)
}
}
}
}
}
}
models <- unique(models)
return(models)
}
## Function to put it all together:
allom_gams <- function(data1, y, x, by, rank = "AIC") {
## y, x and by must be input as concatenated string variables
x <- c("logDBHcm", "logHTcm", "logPVcm3")
by <- c("SPP4", "STANDAGE", "SPPSTAGE")
by <- append(by, " ")
op <- c(" + ", " * ")
## call function to build models for case with one explanatory variable
model_terms <- build_model_terms(x = x, by = by, op = op)
int_terms <- model_terms[[1]]
s_terms <- model_terms[[2]]
models_1 <- build_models_1(y = y, x = x, by = by, int_terms = int_terms, s_terms = s_terms)
models_2 <- build_models_2(y = y, x = x, by = by, int_terms = int_terms, s_terms = s_terms)
models <- c(models_1, models_2)
## fit gam models:
for (i in 1:length(models)) {
assign(paste0("model", i), gam(models[[i]], data = data1))
}
model_list <- character()
for (i in 1:length(models)) {
if (i < max(length(models))) {
model_list <- paste0(model_list, "model", i, ", ")
}
else {
model_list <- paste0(model_list, "model", i)
}
}
results <- eval(parse(text = paste("model.sel(", model_list, ", rank =", rank, ")")))
return(results)
}
## build, run and perform model selection on models
model_sel_gam_AIC <- allom_gams(data1 = allom_data_trunc, y = "logABOVEg", x = c("logDBHcm", "logHTcm", "logPVcm3"), by = c("SPP4", "STANDAGE", "SPPSTAGE"), rank = "AIC")
## Check model selection on dataset without poorly represented species to make sure results don't change:
# model_sel_gam_AICtruncSPP4 <- allom_gams(data1 = allom_data_truncSPP4, y = "logABOVEg", x = c("logDBHcm", "logHTcm", "logPVcm3"), by = c("SPP4", "STANDAGE", "SPPSTAGE"), rank = "AIC")
#
# ## Check model selection on dataset without poorly represented species/standage combos:
# model_sel_gam_AICtruncSPPSTAGE <- allom_gams(data1 = allom_data_truncSPPSTAGE, y = "logABOVEg", x = c("logDBHcm", "logHTcm", "logPVcm3"), by = c("SPP4", "STANDAGE", "SPPSTAGE"), rank = "AIC")
#
# ## Check model selection when STANDAGE treated as numeric instead of factor:
# model_sel_gam_AICtruncSTANDAGENUM <- allom_gams(data1 = allom_data_truncSTANDAGENUM, y = "logABOVEg", x = c("logDBHcm", "logHTcm", "logPVcm3"), by = c("SPP4", "STANDAGE", "SPPSTAGE"), rank = "AIC")Here are the results when using AIC. The output becomes quite cumbersome when linear and gam models are combined because of the number of possible terms between the two model families. The four models with the lowest AIC score are gam models, the fifth being linear.
Here they are listed out in a more readable format:
1) logABOVEg ~ s(logDBHcm, by = SPPSTAGE) + s(logHTcm, by = SPPSTAGE) + SPPSTAGE
2) logABOVEg ~ s(logDBHcm, by = SPPSTAGE) + s(logHTcm, by = SPP4) + SPPSTAGE
3) logABOVEg ~ s(logDBHcm, by = SPPSTAGE) + s(logHTcm, by = STANDAGE) + SPPSTAGE
4) logABOVEg ~ s(logPVcm3, by = SPPSTAGE) + SPPSTAGE
5) logABOVEg ~ logDBHcm + logHTcm + SPP4 + STANDAGE + logDBHcm:SPP4 + logDBH:STANDAGE +
logHT:SPP4 + logHT:STANDAGE + SPP4:STANDAGE + logHT:SPP4:STANDAGE
model_sel_results_AIC <- merge(model_sel_linear_AIC, model_sel_gam_AIC)
head(model_sel_results_AIC)## Model selection table
## (Int) SPP4 STA s(lDB,4,SPPS,1) SPPS s(lHT,4,SPP4,1)
## model164.y 7.949 +
## model162.y 7.867 + +
## model159.y 8.130 + + +
## model116.y 7.949 + +
## model163.y 7.961 + +
## model161.y 7.969 +
## s(lHT,4,SPPS,1) s(lDB,4,SPPS) s(lHT,4,SPP4) s(lHT,4,STA)
## model164.y +
## model162.y + +
## model159.y +
## model116.y
## model163.y + +
## model161.y + +
## s(lHT,4,SPPS) df logLik AIC delta weight
## model164.y + 125 145.774 -41.0 0.00 0.681
## model162.y 103 122.893 -38.2 2.85 0.163
## model159.y 107 126.390 -37.4 3.62 0.111
## model116.y + 123 141.314 -35.1 5.94 0.035
## model163.y 93 109.201 -32.0 9.03 0.007
## model161.y 129 144.594 -29.4 11.61 0.002
## Models ranked by AIC(x)allom_model <- gam(logABOVEg ~ s(logDBHcm, by = SPPSTAGE, k = 4) + s(logHTcm, by = SPPSTAGE, k = 4) + SPP4 * STANDAGE, data = allom_data_trunc, qr=FALSE)
summary(allom_model)##
## Family: gaussian
## Link function: identity
##
## Formula:
## logABOVEg ~ s(logDBHcm, by = SPPSTAGE, k = 4) + s(logHTcm, by = SPPSTAGE,
## k = 4) + SPP4 * STANDAGE
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.91407 0.55737 14.199 <2e-16 ***
## SPP4ACRU 1.06915 0.64682 1.653 0.1000 .
## SPP4ACSA 1.28603 0.63987 2.010 0.0459 *
## SPP4BEAL 1.20918 0.65695 1.841 0.0672 .
## SPP4BEPA 1.22375 0.99314 1.232 0.2194
## SPP4FAGR 1.71842 0.81252 2.115 0.0357 *
## SPP4PIRU -0.02006 0.13693 -0.146 0.8837
## SPP4PRPE 0.90884 0.67657 1.343 0.1808
## STANDAGE16 -5.71597 2.22182 -2.573 0.0108 *
## STANDAGE17 0.03533 0.61209 0.058 0.9540
## STANDAGE23 -1.61731 5.67910 -0.285 0.7761
## STANDAGE26 0.96506 0.55979 1.724 0.0863 .
## STANDAGE29 0.97583 4.09147 0.239 0.8117
## STANDAGE98 0.97960 0.56124 1.745 0.0825 .
## SPP4ACRU:STANDAGE16 5.60162 2.25370 2.486 0.0138 *
## SPP4ACSA:STANDAGE16 4.66599 2.05804 2.267 0.0245 *
## SPP4BEAL:STANDAGE16 4.25476 5.54187 0.768 0.4436
## SPP4BEPA:STANDAGE16 4.61852 2.49501 1.851 0.0657 .
## SPP4FAGR:STANDAGE16 -29.79298 12.04542 -2.473 0.0143 *
## SPP4PIRU:STANDAGE16 0.00000 0.00000 NA NA
## SPP4PRPE:STANDAGE16 4.93611 2.26701 2.177 0.0307 *
## SPP4ACRU:STANDAGE17 0.00000 0.00000 NA NA
## SPP4ACSA:STANDAGE17 4.35475 3.10772 1.401 0.1628
## SPP4BEAL:STANDAGE17 -0.86693 0.78604 -1.103 0.2714
## SPP4BEPA:STANDAGE17 -0.05194 1.02853 -0.050 0.9598
## SPP4FAGR:STANDAGE17 8.50206 11.31072 0.752 0.4532
## SPP4PIRU:STANDAGE17 0.00000 0.00000 NA NA
## SPP4PRPE:STANDAGE17 0.03935 0.72418 0.054 0.9567
## SPP4ACRU:STANDAGE23 0.00000 0.00000 NA NA
## SPP4ACSA:STANDAGE23 0.00000 0.00000 NA NA
## SPP4BEAL:STANDAGE23 2.25911 5.78284 0.391 0.6965
## SPP4BEPA:STANDAGE23 1.49034 5.73253 0.260 0.7952
## SPP4FAGR:STANDAGE23 0.00000 0.00000 NA NA
## SPP4PIRU:STANDAGE23 0.00000 0.00000 NA NA
## SPP4PRPE:STANDAGE23 -5.36675 16.98443 -0.316 0.7524
## SPP4ACRU:STANDAGE26 -1.34815 0.99896 -1.350 0.1788
## SPP4ACSA:STANDAGE26 -1.01868 0.69343 -1.469 0.1435
## SPP4BEAL:STANDAGE26 -0.93070 0.69367 -1.342 0.1813
## SPP4BEPA:STANDAGE26 -1.31017 1.01748 -1.288 0.1994
## SPP4FAGR:STANDAGE26 -1.33278 0.82017 -1.625 0.1058
## SPP4PIRU:STANDAGE26 0.00000 0.00000 NA NA
## SPP4PRPE:STANDAGE26 -0.82925 0.73802 -1.124 0.2626
## SPP4ACRU:STANDAGE29 -0.95117 4.11280 -0.231 0.8174
## SPP4ACSA:STANDAGE29 -0.24322 4.20843 -0.058 0.9540
## SPP4BEAL:STANDAGE29 -0.68423 4.11907 -0.166 0.8682
## SPP4BEPA:STANDAGE29 4.95649 24.49107 0.202 0.8398
## SPP4FAGR:STANDAGE29 -1.45498 4.12750 -0.353 0.7248
## SPP4PIRU:STANDAGE29 0.00000 0.00000 NA NA
## SPP4PRPE:STANDAGE29 -0.64707 4.15499 -0.156 0.8764
## SPP4ACRU:STANDAGE98 0.00000 0.00000 NA NA
## SPP4ACSA:STANDAGE98 -1.15923 0.64601 -1.794 0.0743 .
## SPP4BEAL:STANDAGE98 -1.07396 0.66310 -1.620 0.1070
## SPP4BEPA:STANDAGE98 0.00000 0.00000 NA NA
## SPP4FAGR:STANDAGE98 -1.49312 0.81746 -1.827 0.0693 .
## SPP4PIRU:STANDAGE98 -0.02006 0.13693 -0.146 0.8837
## SPP4PRPE:STANDAGE98 0.00000 0.00000 NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(logDBHcm):SPPSTAGEACPE 17 1.0000 1.0000 12.602 0.000482 ***
## s(logDBHcm):SPPSTAGEACPE 26 1.0000 1.0000 28.241 2.81e-07 ***
## s(logDBHcm):SPPSTAGEACPE 98 1.0000 1.0000 28.734 2.25e-07 ***
## s(logDBHcm):SPPSTAGEACRU 14 1.0000 1.0000 14.472 0.000190 ***
## s(logDBHcm):SPPSTAGEACRU 16 1.0000 1.0000 6.096 0.014402 *
## s(logDBHcm):SPPSTAGEACRU 26 1.0000 1.0000 2.735 0.099759 .
## s(logDBHcm):SPPSTAGEACRU 29 1.0000 1.0000 4.961 0.027058 *
## s(logDBHcm):SPPSTAGEACSA 14 1.0000 1.0000 18.927 2.17e-05 ***
## s(logDBHcm):SPPSTAGEACSA 16 0.4606 0.4606 0.025 0.915474
## s(logDBHcm):SPPSTAGEACSA 17 1.0000 1.0000 10.772 0.001219 **
## s(logDBHcm):SPPSTAGEACSA 26 1.0000 1.0000 17.680 3.95e-05 ***
## s(logDBHcm):SPPSTAGEACSA 29 1.0000 1.0000 1.589 0.208941
## s(logDBHcm):SPPSTAGEACSA 98 1.4033 1.6758 155.808 < 2e-16 ***
## s(logDBHcm):SPPSTAGEBEAL 14 1.0000 1.0000 8.167 0.004727 **
## s(logDBHcm):SPPSTAGEBEAL 16 1.0000 1.0000 0.038 0.845011
## s(logDBHcm):SPPSTAGEBEAL 17 1.0000 1.0000 14.943 0.000150 ***
## s(logDBHcm):SPPSTAGEBEAL 23 1.0000 1.0000 1.524 0.218468
## s(logDBHcm):SPPSTAGEBEAL 26 1.0000 1.0000 13.976 0.000242 ***
## s(logDBHcm):SPPSTAGEBEAL 29 1.0000 1.0000 4.091 0.044488 *
## s(logDBHcm):SPPSTAGEBEAL 98 1.4009 1.6002 197.387 < 2e-16 ***
## s(logDBHcm):SPPSTAGEBEPA 14 1.0000 1.0000 2.361 0.126018
## s(logDBHcm):SPPSTAGEBEPA 16 1.0000 1.0000 0.255 0.614217
## s(logDBHcm):SPPSTAGEBEPA 17 1.0000 1.0000 50.511 1.67e-11 ***
## s(logDBHcm):SPPSTAGEBEPA 23 1.0000 1.0000 26.291 6.86e-07 ***
## s(logDBHcm):SPPSTAGEBEPA 26 1.0000 1.0000 51.077 1.31e-11 ***
## s(logDBHcm):SPPSTAGEBEPA 29 1.0000 1.0000 0.078 0.780913
## s(logDBHcm):SPPSTAGEFAGR 14 1.0000 1.0000 6.648 0.010661 *
## s(logDBHcm):SPPSTAGEFAGR 16 1.0000 1.0000 5.996 0.015222 *
## s(logDBHcm):SPPSTAGEFAGR 17 1.6803 1.8858 23.388 2.36e-07 ***
## s(logDBHcm):SPPSTAGEFAGR 26 1.0000 1.0000 39.876 1.59e-09 ***
## s(logDBHcm):SPPSTAGEFAGR 29 0.3243 0.3243 61.359 1.36e-05 ***
## s(logDBHcm):SPPSTAGEFAGR 98 1.0000 1.0000 361.507 < 2e-16 ***
## s(logDBHcm):SPPSTAGEPIRU 98 1.0000 1.0000 36.534 6.87e-09 ***
## s(logDBHcm):SPPSTAGEPRPE 14 1.0000 1.0000 5.816 0.016800 *
## s(logDBHcm):SPPSTAGEPRPE 16 1.0000 1.0000 3.494 0.063084 .
## s(logDBHcm):SPPSTAGEPRPE 17 1.0000 1.0000 23.047 3.09e-06 ***
## s(logDBHcm):SPPSTAGEPRPE 23 1.5500 1.7942 2.852 0.126479
## s(logDBHcm):SPPSTAGEPRPE 26 1.0000 1.0000 23.426 2.59e-06 ***
## s(logDBHcm):SPPSTAGEPRPE 29 1.0000 1.0000 3.862 0.050810 .
## s(logHTcm):SPPSTAGEACPE 17 1.0000 1.0000 0.551 0.458696
## s(logHTcm):SPPSTAGEACPE 26 1.0000 1.0000 0.042 0.838394
## s(logHTcm):SPPSTAGEACPE 98 1.0000 1.0000 0.019 0.890443
## s(logHTcm):SPPSTAGEACRU 14 1.0000 1.0000 0.426 0.514697
## s(logHTcm):SPPSTAGEACRU 16 1.0000 1.0000 0.128 0.720471
## s(logHTcm):SPPSTAGEACRU 26 1.0000 1.0000 0.437 0.509197
## s(logHTcm):SPPSTAGEACRU 29 1.0000 1.0000 0.000 0.986487
## s(logHTcm):SPPSTAGEACSA 14 1.0000 1.0000 0.006 0.938250
## s(logHTcm):SPPSTAGEACSA 16 0.6146 0.6146 35.906 4.89e-06 ***
## s(logHTcm):SPPSTAGEACSA 17 1.8681 1.9836 3.303 0.052826 .
## s(logHTcm):SPPSTAGEACSA 26 1.0000 1.0000 0.057 0.812109
## s(logHTcm):SPPSTAGEACSA 29 1.0000 1.0000 0.295 0.587747
## s(logHTcm):SPPSTAGEACSA 98 1.0000 1.0000 0.713 0.399588
## s(logHTcm):SPPSTAGEBEAL 14 1.0000 1.0000 1.106 0.294254
## s(logHTcm):SPPSTAGEBEAL 16 1.0000 1.0000 0.101 0.750507
## s(logHTcm):SPPSTAGEBEAL 17 2.8624 2.9821 12.693 4.03e-07 ***
## s(logHTcm):SPPSTAGEBEAL 23 1.0000 1.0000 0.081 0.776357
## s(logHTcm):SPPSTAGEBEAL 26 1.5280 1.8050 1.201 0.399854
## s(logHTcm):SPPSTAGEBEAL 29 1.0000 1.0000 0.097 0.755601
## s(logHTcm):SPPSTAGEBEAL 98 1.5447 1.7934 0.799 0.345325
## s(logHTcm):SPPSTAGEBEPA 14 1.0000 1.0000 0.125 0.724112
## s(logHTcm):SPPSTAGEBEPA 16 1.0000 1.0000 2.010 0.157857
## s(logHTcm):SPPSTAGEBEPA 17 1.1492 1.2780 0.259 0.584128
## s(logHTcm):SPPSTAGEBEPA 23 1.0000 1.0000 0.065 0.799534
## s(logHTcm):SPPSTAGEBEPA 26 1.0000 1.0000 0.492 0.483732
## s(logHTcm):SPPSTAGEBEPA 29 1.0000 1.0000 0.044 0.834086
## s(logHTcm):SPPSTAGEFAGR 14 1.0000 1.0000 0.606 0.437345
## s(logHTcm):SPPSTAGEFAGR 16 1.0000 1.0000 6.111 0.014285 *
## s(logHTcm):SPPSTAGEFAGR 17 1.1686 1.2967 1.387 0.415472
## s(logHTcm):SPPSTAGEFAGR 26 1.0000 1.0000 0.043 0.835621
## s(logHTcm):SPPSTAGEFAGR 29 0.6781 0.6781 62.894 4.79e-10 ***
## s(logHTcm):SPPSTAGEFAGR 98 2.0526 2.2926 17.211 2.04e-07 ***
## s(logHTcm):SPPSTAGEPIRU 98 2.5078 2.8081 2.889 0.024406 *
## s(logHTcm):SPPSTAGEPRPE 14 1.0000 1.0000 0.194 0.660419
## s(logHTcm):SPPSTAGEPRPE 16 1.0000 1.0000 31.753 5.75e-08 ***
## s(logHTcm):SPPSTAGEPRPE 17 1.0000 1.0000 41.740 7.09e-10 ***
## s(logHTcm):SPPSTAGEPRPE 23 1.0061 1.0088 0.065 0.851093
## s(logHTcm):SPPSTAGEPRPE 26 1.5285 1.7808 0.313 0.720657
## s(logHTcm):SPPSTAGEPRPE 29 1.0000 1.0000 0.134 0.715177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Rank: 271/290
## R-sq.(adj) = 0.994 Deviance explained = 99.6%
## GCV = 0.063204 Scale est. = 0.038352 n = 316Plotting the model:
plot(allom_model, scale=0)
## doesn't seem to be any discernible pattern, however it is hard to determine with 42 plots
## Plotting predicted values, observed values, and confidence interval for each combination of spp and standage:Write function to plot predictions vs observed for all combinations of SPP and STANDAGE:
plot_pred <- function(data, model, pred.length = 50, plotcov = c("logDBHcm, logHTcm")) {
## Build predicted values into dataframe:
numspp <- length(levels(data$SPP4))
numstage <- length(levels(data$STANDAGE))
## write in STANDAGE levels for each SPP
stagevec <- rep(levels(data$STANDAGE), numspp)
## write in SPP
sppvec <- character()
for (i in 1:numspp) {
newspp <- rep(levels(data$SPP4)[i], numstage)
sppvec <- c(sppvec, newspp)
}
## combine into dataframe
fac_data <- data.frame(sppvec, stagevec)
## get x scale for each STANDAGE and covariate
STAGE <- levels(data$STANDAGE)
minDBH <- numeric(length = length(levels(data$STANDAGE)))
maxDBH <- numeric(length = length(levels(data$STANDAGE)))
minHT <- numeric(length = length(levels(data$STANDAGE)))
maxHT <- numeric(length = length(levels(data$STANDAGE)))
ranges <- data.frame(STAGE, minDBH, maxDBH, minHT, maxHT)
for (i in 1:length(levels(data$STANDAGE))) {
newsub <- subset(data, STANDAGE == levels(data$STANDAGE)[i])
ranges$minDBH[i] <- min(newsub$logDBHcm)
ranges$maxDBH[i] <- max(newsub$logDBHcm)
ranges$minHT[i] <- min(newsub$logHTcm)
ranges$maxHT[i] <- max(newsub$logHTcm)
}
## write in regularly spaces covariate values to predict upon
logDBHcm <- numeric()
logHTcm <- numeric()
for (i in 1:nrow(fac_data)) {
range <- subset(ranges, STAGE == stagevec[i])
logDBHcm <- c(logDBHcm, seq(range$minDBH, range$maxDBH, length.out = pred.length))
logHTcm <- c(logHTcm, seq(range$minHT, range$maxHT, length.out = pred.length))
}
SPP4 <- rep(fac_data$sppvec, each = pred.length)
STANDAGE <- rep(fac_data$stagevec, each = pred.length)
new_data <- data.frame(SPP4, STANDAGE, logDBHcm, logHTcm)
## Create SPPSTAGE variable in new data for prediction
for (i in 1:length(new_data$SPP4)) {
new_data$SPPSTAGE[i] <- paste(new_data$SPP4[i], new_data$STANDAGE[i])
}
new_data$SPPSTAGE <- as.factor(new_data$SPPSTAGE)
## Eliminate covariate combinations in predictive data for which there was no observed data
del_vec <- character()
for (i in 1:nrow(new_data)) {
if (!any(grepl(new_data$SPPSTAGE[i], levels(data$SPPSTAGE)))) {
new <- rownames(new_data)[i]
del_vec <- c(del_vec, new)
}
else {}
}
del_vec <- as.numeric(del_vec[!is.na(del_vec)])
new_data <- new_data[-del_vec,]
## write predictions to object
predicted_data <- predict(model, newdata = new_data, se.fit = TRUE)
## write in predictions and standard errors to dataframe
new_data$logABOVEg <- predicted_data$fit
new_data$SE <- predicted_data$se.fit
new_data$CI.low <- new_data$logABOVEg - predicted_data$se.fit
new_data$CI.high <- new_data$logABOVEg + predicted_data$se.fit
plots <- list()
if (grepl("logDBHcm", plotcov) && grepl("logHTcm", plotcov)) {
DBHplot <- ggplot(new_data, aes(logDBHcm, logABOVEg)) +
geom_line(size = 1) +
geom_line(aes(logDBHcm, CI.high), color = "palevioletred1") +
geom_line(aes(logDBHcm, CI.low), color = "palevioletred1") +
geom_jitter(data = data, aes(logDBHcm, logABOVEg), size = 0.75) +
facet_grid(SPP4 ~ STANDAGE, scales = "free")
HTplot <- ggplot(new_data, aes(logHTcm, logABOVEg)) +
geom_line(size = 1) +
geom_line(aes(logHTcm, CI.high), color = "palevioletred1") +
geom_line(aes(logHTcm, CI.low), color = "palevioletred1") +
geom_jitter(data = data, aes(logHTcm, logABOVEg), size = 0.75) +
facet_grid(SPP4 ~ STANDAGE, scales = "free")
print(DBHplot)
print(HTplot)
}
else if (grepl("logHTcm", plotcov) && !grepl("logDBHcm", plotcov)) {
HTplot <- ggplot(new_data, aes(logHTcm, logABOVEg)) +
geom_line(size = 1) +
geom_line(aes(logHTcm, CI.high), color = "palevioletred1") +
geom_line(aes(logHTcm, CI.low), color = "palevioletred1") +
geom_jitter(data = allom_data_trunc, aes(logHTcm, logABOVEg), size = 0.75) +
facet_grid(SPP4 ~ STANDAGE, scales = "free")
print(HTplot)
}
else if (grepl("logDBHcm", plotcov) && !grepl("logHTcm", plotcov)) {
DBHplot <- ggplot(new_data, aes(logDBHcm, logABOVEg)) +
geom_line(size = 1) +
geom_line(aes(logDBHcm, CI.high), color = "palevioletred1") +
geom_line(aes(logDBHcm, CI.low), color = "palevioletred1") +
geom_jitter(data = data, aes(logDBHcm, logABOVEg), size = 0.75) +
facet_grid(SPP4 ~ STANDAGE, scales = "free")
print(DBHplot)
}
else {
print("error: plot cov must contain either logHTcm, logDBHcm or both")
}
}
plot_pred(data = allom_data_trunc, model = allom_model)
anova(allom_model)##
## Family: gaussian
## Link function: identity
##
## Formula:
## logABOVEg ~ s(logDBHcm, by = SPPSTAGE, k = 4) + s(logHTcm, by = SPPSTAGE,
## k = 4) + SPP4 * STANDAGE
##
## Parametric Terms:
## df F p-value
## SPP4 7 0.992 0.438
## STANDAGE 6 2.578 0.020
## SPP4:STANDAGE 30 4.586 3.36e-11
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(logDBHcm):SPPSTAGEACPE 17 1.0000 1.0000 12.602 0.000482
## s(logDBHcm):SPPSTAGEACPE 26 1.0000 1.0000 28.241 2.81e-07
## s(logDBHcm):SPPSTAGEACPE 98 1.0000 1.0000 28.734 2.25e-07
## s(logDBHcm):SPPSTAGEACRU 14 1.0000 1.0000 14.472 0.000190
## s(logDBHcm):SPPSTAGEACRU 16 1.0000 1.0000 6.096 0.014402
## s(logDBHcm):SPPSTAGEACRU 26 1.0000 1.0000 2.735 0.099759
## s(logDBHcm):SPPSTAGEACRU 29 1.0000 1.0000 4.961 0.027058
## s(logDBHcm):SPPSTAGEACSA 14 1.0000 1.0000 18.927 2.17e-05
## s(logDBHcm):SPPSTAGEACSA 16 0.4606 0.4606 0.025 0.915474
## s(logDBHcm):SPPSTAGEACSA 17 1.0000 1.0000 10.772 0.001219
## s(logDBHcm):SPPSTAGEACSA 26 1.0000 1.0000 17.680 3.95e-05
## s(logDBHcm):SPPSTAGEACSA 29 1.0000 1.0000 1.589 0.208941
## s(logDBHcm):SPPSTAGEACSA 98 1.4033 1.6758 155.808 < 2e-16
## s(logDBHcm):SPPSTAGEBEAL 14 1.0000 1.0000 8.167 0.004727
## s(logDBHcm):SPPSTAGEBEAL 16 1.0000 1.0000 0.038 0.845011
## s(logDBHcm):SPPSTAGEBEAL 17 1.0000 1.0000 14.943 0.000150
## s(logDBHcm):SPPSTAGEBEAL 23 1.0000 1.0000 1.524 0.218468
## s(logDBHcm):SPPSTAGEBEAL 26 1.0000 1.0000 13.976 0.000242
## s(logDBHcm):SPPSTAGEBEAL 29 1.0000 1.0000 4.091 0.044488
## s(logDBHcm):SPPSTAGEBEAL 98 1.4009 1.6002 197.387 < 2e-16
## s(logDBHcm):SPPSTAGEBEPA 14 1.0000 1.0000 2.361 0.126018
## s(logDBHcm):SPPSTAGEBEPA 16 1.0000 1.0000 0.255 0.614217
## s(logDBHcm):SPPSTAGEBEPA 17 1.0000 1.0000 50.511 1.67e-11
## s(logDBHcm):SPPSTAGEBEPA 23 1.0000 1.0000 26.291 6.86e-07
## s(logDBHcm):SPPSTAGEBEPA 26 1.0000 1.0000 51.077 1.31e-11
## s(logDBHcm):SPPSTAGEBEPA 29 1.0000 1.0000 0.078 0.780913
## s(logDBHcm):SPPSTAGEFAGR 14 1.0000 1.0000 6.648 0.010661
## s(logDBHcm):SPPSTAGEFAGR 16 1.0000 1.0000 5.996 0.015222
## s(logDBHcm):SPPSTAGEFAGR 17 1.6803 1.8858 23.388 2.36e-07
## s(logDBHcm):SPPSTAGEFAGR 26 1.0000 1.0000 39.876 1.59e-09
## s(logDBHcm):SPPSTAGEFAGR 29 0.3243 0.3243 61.359 1.36e-05
## s(logDBHcm):SPPSTAGEFAGR 98 1.0000 1.0000 361.507 < 2e-16
## s(logDBHcm):SPPSTAGEPIRU 98 1.0000 1.0000 36.534 6.87e-09
## s(logDBHcm):SPPSTAGEPRPE 14 1.0000 1.0000 5.816 0.016800
## s(logDBHcm):SPPSTAGEPRPE 16 1.0000 1.0000 3.494 0.063084
## s(logDBHcm):SPPSTAGEPRPE 17 1.0000 1.0000 23.047 3.09e-06
## s(logDBHcm):SPPSTAGEPRPE 23 1.5500 1.7942 2.852 0.126479
## s(logDBHcm):SPPSTAGEPRPE 26 1.0000 1.0000 23.426 2.59e-06
## s(logDBHcm):SPPSTAGEPRPE 29 1.0000 1.0000 3.862 0.050810
## s(logHTcm):SPPSTAGEACPE 17 1.0000 1.0000 0.551 0.458696
## s(logHTcm):SPPSTAGEACPE 26 1.0000 1.0000 0.042 0.838394
## s(logHTcm):SPPSTAGEACPE 98 1.0000 1.0000 0.019 0.890443
## s(logHTcm):SPPSTAGEACRU 14 1.0000 1.0000 0.426 0.514697
## s(logHTcm):SPPSTAGEACRU 16 1.0000 1.0000 0.128 0.720471
## s(logHTcm):SPPSTAGEACRU 26 1.0000 1.0000 0.437 0.509197
## s(logHTcm):SPPSTAGEACRU 29 1.0000 1.0000 0.000 0.986487
## s(logHTcm):SPPSTAGEACSA 14 1.0000 1.0000 0.006 0.938250
## s(logHTcm):SPPSTAGEACSA 16 0.6146 0.6146 35.906 4.89e-06
## s(logHTcm):SPPSTAGEACSA 17 1.8681 1.9836 3.303 0.052826
## s(logHTcm):SPPSTAGEACSA 26 1.0000 1.0000 0.057 0.812109
## s(logHTcm):SPPSTAGEACSA 29 1.0000 1.0000 0.295 0.587747
## s(logHTcm):SPPSTAGEACSA 98 1.0000 1.0000 0.713 0.399588
## s(logHTcm):SPPSTAGEBEAL 14 1.0000 1.0000 1.106 0.294254
## s(logHTcm):SPPSTAGEBEAL 16 1.0000 1.0000 0.101 0.750507
## s(logHTcm):SPPSTAGEBEAL 17 2.8624 2.9821 12.693 4.03e-07
## s(logHTcm):SPPSTAGEBEAL 23 1.0000 1.0000 0.081 0.776357
## s(logHTcm):SPPSTAGEBEAL 26 1.5280 1.8050 1.201 0.399854
## s(logHTcm):SPPSTAGEBEAL 29 1.0000 1.0000 0.097 0.755601
## s(logHTcm):SPPSTAGEBEAL 98 1.5447 1.7934 0.799 0.345325
## s(logHTcm):SPPSTAGEBEPA 14 1.0000 1.0000 0.125 0.724112
## s(logHTcm):SPPSTAGEBEPA 16 1.0000 1.0000 2.010 0.157857
## s(logHTcm):SPPSTAGEBEPA 17 1.1492 1.2780 0.259 0.584128
## s(logHTcm):SPPSTAGEBEPA 23 1.0000 1.0000 0.065 0.799534
## s(logHTcm):SPPSTAGEBEPA 26 1.0000 1.0000 0.492 0.483732
## s(logHTcm):SPPSTAGEBEPA 29 1.0000 1.0000 0.044 0.834086
## s(logHTcm):SPPSTAGEFAGR 14 1.0000 1.0000 0.606 0.437345
## s(logHTcm):SPPSTAGEFAGR 16 1.0000 1.0000 6.111 0.014285
## s(logHTcm):SPPSTAGEFAGR 17 1.1686 1.2967 1.387 0.415472
## s(logHTcm):SPPSTAGEFAGR 26 1.0000 1.0000 0.043 0.835621
## s(logHTcm):SPPSTAGEFAGR 29 0.6781 0.6781 62.894 4.79e-10
## s(logHTcm):SPPSTAGEFAGR 98 2.0526 2.2926 17.211 2.04e-07
## s(logHTcm):SPPSTAGEPIRU 98 2.5078 2.8081 2.889 0.024406
## s(logHTcm):SPPSTAGEPRPE 14 1.0000 1.0000 0.194 0.660419
## s(logHTcm):SPPSTAGEPRPE 16 1.0000 1.0000 31.753 5.75e-08
## s(logHTcm):SPPSTAGEPRPE 17 1.0000 1.0000 41.740 7.09e-10
## s(logHTcm):SPPSTAGEPRPE 23 1.0061 1.0088 0.065 0.851093
## s(logHTcm):SPPSTAGEPRPE 26 1.5285 1.7808 0.313 0.720657
## s(logHTcm):SPPSTAGEPRPE 29 1.0000 1.0000 0.134 0.715177## looks like the SPP4 parametric terms arent significant. STANDAGE and interactions are both significant.
allom_p.coeff <- data.frame(summary(allom_model)$p.coeff)
coefnames <- as.factor(attributes(summary(allom_model)$p.coeff)$names)
coefvals <- summary(allom_model)$p.coeff
attributes(coefvals) <- NULL
allom_coef <- data.frame(coefnames, coefvals)
allom_coefSTANDAGE <- allom_coef[9:14, ]
allom_coefSPP4 <- allom_coef[2:8, ]
## there appears to be some sort of pattern here
ggplot(allom_coefSTANDAGE, aes(x = coefnames, y = coefvals)) +
geom_point() 
ggplot(allom_coefSPP4, aes(x = coefnames, y = coefvals)) +
geom_point()
## x axis is illegible, but could be useful for seeing larger trends
ggplot(allom_coef, aes(x = coefnames, y = coefvals)) +
geom_point()
with STANDAGE as class numeric (IGNORE FOR NOW)
# allom_p.coeff <- data.frame(summary(allom_model)$p.coeff)
#
# coefnames <- as.factor(attributes(summary(allom_model)$p.coeff)$names)
# coefvals <- summary(allom_model)$p.coeff
# attributes(coefvals) <- NULL
#
# allom_coef <- data.frame(coefnames, coefvals)
#
# allom_coefSPP4 <- allom_coef[2:8, ]
# allom_coefSTANDAGE <- allom_coef[9:14, ]
# allom_coefSTANDAGE$coefnames <- as.character(allom_coefSTANDAGE$coefnames)
#
# for (i in 1:nrow(allom_coefSTANDAGE)) {
# allom_coefSTANDAGE[i, 1] <- gsub("STANDAGE", "", allom_coefSTANDAGE[i, 1])
# }
#
# allom_coefSTANDAGE$coefnames <- as.numeric(allom_coefSTANDAGE$coefnames)
# allom_coefSTANDAGE
#
# ggplot(allom_coefSTANDAGE, aes(x = coefnames, y = coefvals)) +
# geom_point()
#
# ggplot(allom_coefSPP4, aes(x = coefnames, y = coefvals)) +
# geom_point()
#
# ## x axis is illegible, but could be useful for seeing larger trends
# ggplot(allom_coef, aes(x = coefnames, y = coefvals)) +
# geom_point()